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Abstract, 

The  main  objective  of  this  study  is  to  develop  an  efficient  multiscale  coarse  grid  method  which  can  ' 
be  used  as  a  competitive  algorithm  in  studying  composite  materials  and  flow  transport  in  strongly 
heterogeneous  porous  media.  On  one  hand,  we  have  explored  the  possibility  of  using  adaptive 
mesh  to  reduce  the  modeling  error  introduced  by  the  traditional  moment  average  technique.  On 
the  other  hand,  we  found  that  in  the  case  of  high  aqeet*  ratio  permeability  tensor,  the  modeling 
error  in  ignoring  high  order  moments  (3rd  order  or  higher)  could  be  very  large.  To  overcome 
this  difficulty,  we  have  investigated  an  alternative  approach  that  uses  two-scale  homogenization 
analysis  to  derive  a  coarse  grid  model  in  a  systematic  way.  Finally,  we  have  made  some  progress  in 
developing  numerical  methods  to  solve  multiscale  nonlinear  stochastic  partial  differential  equations  , 
by  using  Wiener-Chaos  expansions.  These  methods  will  reduce  the  problem  of  solving  stochastic 
PDEs  to  solving  a  set  of  deterministic  PDEs,  This  numerical  method  can  be  combined  with  our 
multiscale  computational  method,  and  can  be  used  to  compute  accurately  high  order  statistical 
quantites  more  efficiently  than  the  traditional  Monte-Carlo  method. 
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I.  Technical  Accomplishments. 

(1)  The  over-sampling  technique  and  its  error  analysis. 

One  of  the  main  difficulties  in  developing  the  multiscale  computational  methods  is  how  to  ef¬ 
fectively  localize  the  equation  governing  the  small  scales.  If  we  use  artificial  boundary  conditions 
from  the  large  scales,  such  as  linear  boundary  conditions,  it  will  create  mismatch  among  the  small 
scales  across  the  boundaries  of  coarse  grid  elements.  Our  analysis  [6,  7]  indicates  that  the  error  . 
appears  as  a  resonance  between  the  small  physical  scales  of  the  medium  and  the  mesh  size.  Indeed 
the  error  is  given  by  the  ratio  between  the  two  scales;  it  increases  as  the  size  of  grid  blocks  gets 
close  to  the  small  physical  scales.  We  also  show  that  the  effect  of  different  boundary  conditions  lies 
in  a  narrow  region  near  the  boundaries  of  grid  blocks  and  it  contributes  to  part  of  the  resonance 
error.  In  [6],  we  introduced  an  over-sampling  technique  to  remove  the  boundary  layer  effect.  This 
is  a  very  effective  method  that  can  be  implemented  easily  for  problems  with  many  or  continuous 
scales.  This  idea  has  recently  been  modified  for  the  finite  volume  method  and  used  in  the  new  flow 
simulator  code  of  Chevron. 

A  consequence  of  the  over-sampling  method  is  that  the  resulting  multiscale  finite  element  method 
is  no  longer  conforming.  With  Efendiev  and  Wu  [5],  we  give  a  detailed  analysis  of  the  non- 
conforming  error.  Our  analysis  also  reveals  a  new  cell  resonance  error  which  is  caused  by  the 
mismatch  between  the  mesh  size  and  the  wavelength  of  the  small  scale.  We  show  that  the  cell 
resonance  error  is  of  lower  order  and  is  difficult  to  observe  in  practice. 

(2) .  A  mixed  multiscale  finite  element  method. 

Recently,  in  collaboration  with  Dr.  Z,  Chen  [4],  we  have  developed  a  mixed  multiscale  finite 
element  method.  The  advantage  of  a  mixed  finite  element  formulation  is  that  the  velocity  flux 
is  locally  conserved  across  element  boundaries.  This  is  a  very  important  property  in  many  flow 


simulations  for  large  times.  The  violation  of  this  local  conservation  property  will  lead  to  leakage 
of  velocity  flux.  This  will  deteriorate  the  accuracy  of  the  numerical  solution  for  long  time  compu¬ 
tations.  This  is  the  reason  why  mixed  finite  element  methods  are  attractive  for  porous  medium 
simulations.  Our  computational  results  have  demonstrated  convincingly  that  the  mixed  multiscale 
finite  element  method  gives  more  accurate  results  for  long  time  computations  than  the  displacement 
multiscale  finite  element  method. 

(3) .  A  PDE-based  adaptive  mesh  generator. 

In  collaboration  with  Dr,  Ceniceros  [3],  we  has  developed  a  new  approach  to  generate  a  dynam¬ 
ically  adaptive  mesh  for  flows  with  singular  layered  structures.  The  use  of  well-resolved  uniform 
meshes  for  this  type  of  flows  becomes  prohibitively  expensive.  How  to  construct  an  effective  moving 
meshes  that  can  follow  these  singular  layered  structures  dynamically  is  a  very  challenging  problem. 
By  solving  a  set  of  nonlinear  elliptic  equations,  the  mesh  map  is  generated  in  a  single  step.  The 
resulting  adaptive  method  has  the  following  attractive  properties:  (1)  It  is  fast,  with  an  optimal 
operation  count  O(N),  N  is  the  number  of  grid  points;  (2)  It  can  capture  the  dynamical  evolution  of 
singular  structures  without  introducing  artificial  numerical  diffusion;  (3)  It  is  highly  accurate  and 
stable.  (4)  It  can  be  coupled  easily  with  any  existing  PDE  solver.  We  apply  this  dynamical  moving 
mesh  technique  to  study  the  evolution  of  an  unstably  stratified  flow  with  three  constant  regions  of 
densities  connected  by  two  thin  layers.  This  is  an  extremely  challenging  problem.  The  flow  develops 
very  complex  structure  dynamically  due  to  the  Ravleigh-Taylor  instability.  The  mesh  effectively 
follows  the  flow  throughout  its  complex  evolution.  The  method  is  very  stable  and  efficient  during 
the  entire  computation, 

(4) .  Combining  multiscale  modeling  with  mesh  adaptivity 

We  have  investigated  the  possibility  of  combining  the  multiscale  finite  element  method  with  our 
PDE-based  adaptive  mesh  generator.  The  adaptive  mesh  will  be  used  to  capture  the  dominant  flow 
features,  and  the  multiscale  method  will  be  used  to  capture  the  small  scale  effect  within  each  coarse 
grid  element  by  using  moment  average  techniques.  This  will  provide  an  effective  coarse  grid  method 
for  two-phase  flows  in  heterogeneous  media.  My  Ph.D,  student,  Andy  Westhead,  is  now  working 
at  Exxon-Mobil  as  a  interm  for  the  second  time  to  test  this  idea  for  realistic  flow  simulations.  We 
are  hopeful  that  this  will  lead  to  an  improved  accuracy  in  the  upscaling  model. 

(5) .  Two-scale  analysis  for  incompressible  flows 

A  very  common  phenomenon  in  mechanics,  physics,  chemistry  and  engineering  is  that  processes 
contain  a  wide  range  of  spatial  and  temporal  scales,  which  lead  to  rapidly  varying  structure  in  space 
and  time.  These  include  phase  transitions,  porous  media  flows,  composite  materials,  acoustic  waves 
and  turbulent.  The  analysis  of  flows  with  rapidly  varying  structure  is  a  very  complex  mathematical 
problem.  Direct  numerical  simulations  of  these  problems  require  tremendous  amount  of  computer 
memory  and  CPU  time  which  can  easily  exceed  the  limits  of  modem  computing  resources.  There¬ 
fore,  it  is  desirable  to  develop  numerical  methods  which  capture  the  effect  of  small  scales  on  large 
scales  using  a  relatively  coarse  grid. 

Multiscale  analysis  for  periodic  structures  (homogenization)  has  emerged  as  one  of  the  successful 
techniques  for  these  difficult  problems  with  rapidly  varying  structure.  The  simpliest  multiscale 
analysis  involves  two  scales.  The  first  scale  describes  macroscopic  quantities  or  slow  variables.  The 
second  scale  describes  microscopic  quantities  or  fast  variables.  Homogenization  theory  studies  the 
relation  between  microscopic  and  macroscopic  scales  and  provides  effective  equationes  for  macro¬ 
scopic  quantities.  It  has  been  very  successful  in  solving  elliptic  and  parabolical  problems  with 
rapidly  oscillating  coefficients.  But  ther  has  been  only  limited  success  in  applying  homogenization 
theory  to  hyperbolic  problems. 
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Analysis  of  propagation  of  oscillations  in  incompressible  Euler  and  Navier-Stokes  equations  has 
proved  to  be  an  extremely  ehanllenging  problem.  The  pioneering  work  in  this  area  was  done  by 
Mclaughlin,  Papanicolaou  and  Pironneau  in  their  1985  paper  [12].  By  using  multiple  scale  analysis, 
they  tried  to  obtain  the  homogenized  equation  for  the  incompressible  Euler  equation  under  the 
assumption  that  the  highly  oscillating  part  of  the  velocity  field  is  transported  only  by  the  mean 
flow.  This  crucial  assumption  simplifies  their  multiscale  analysis.  However,  this  assumption  is  not 
very  physical.  In  fact,  there  is  no  uniqueness  of  the  cell  problem  they  formulated.  The  existence  of 
the  cell  problem  is  also  open.  Although  a  number  of  researchers  have  tried  to  resolve  this  difficulty, 
it  is  still  one  of  the  outstanding  open  questions  in  applied  mathematics. 

Recently,  we  have  made  some  progress  in  deriving  a  homogenized  equation  for  incompressible 
Euler  equations  in  two  and  three  space  dimensions.  The  key  idea  is  to  perform  multiscale  expansions 
in  the  Lagrangian  variable.  The  incompressible  Euler  equation  is  reduced  to  a  nonlinear  coupled 
system  of  the  transport  equation,  which  characterizes  the  flow  map,  and  the  elliptic  problem  for 
the  stream  function.  Our  multiscale  analysis  is  carried  out  using  this  framework.  This  nonlinear 
coupled  system  decribes  the  relation  between  large  scales  and  small  scales.  It  provides  a  beautiful 
mathematical  model  to  understand  complicated  phenomenon  of  flows. 

Our  homogenization  theory  for  the  incompressible  Euler  equation  with  multi-scale  initial  data 
shows  that  the  effect  of  small  scale  has  only  a  local  domain  of  influence.  This  enables  us  to 
develop  an  effective  multiscale  numerical  method  based  on  our  asymptotic  analysis  which  captures 
the  correct  large  scale  solution  without  making  closure  assumption.  We  plan  to  generalize  our 
results  to  more  general  multiscale  initial  data  without  restrictive  assumption  on  periodic  structure 
and  scale  separation.  We  will  also  perform  a  number  of  numerical  experiments  to  compare  our 
results  with  the  well-resolved  solution  and  with  other  existing  methods.  Ultimately,  we  want  to 
demonstrate  that  this  two-scale  analysis  can  be  used  to  capture  the  large  time  bebavier  of  the 
macroscopic  solution  for  incompresssible  Navier-Stokes  equations.  If  true,  this  will  provide  an 
effective  numerical  method  for  computing  multiscale  solution  of  incompressible  flows. 

(8).  Numerical  methods  for  solving  stochastic  PDEs 

In  practice,  there  are  many  physical  processes  that  axe  governed  by  stochastic  PDEs.  Examples 
include  wave  propagation  and  diffusion  through  heterogeneous  random  media  and  the  randomly 
forced  Navier-Stokes  equations.  One  commonly  used  method  is  the  Monte  Carlo  method  in  which 
we  simulate  the  problem  realization  by  realization,  and  then  average  the  solutions  over  many 
realizations.  In  each  realization,  one  has  to  use  a  fine  mesh  to  resolve  the  smallest  scales  introduced 
by  the  random  perturbation.  And  one  needs  to  perform  many  realizations  in  order  to  obtain 
accurate  approximations  to  various  statistical  properties.  The  number  of  realizations  could  range 
from  a  few  thousands  to  a  few  hundred  thousands  for  problems  with  large  variance. 

Here  we  are  interested  in  designing  new  type  of  numerical  methods  based  on  some  well-known 
stochastic  decompositions  such  as  Karhunen-Loeve  expansion  [9],  and  Wiener-Chaos  expansion  [2]. 
Wiener-Chaos  expansion  represents  any  function  u(t,  a;,  Xt)  depending  on  a  stochastic  process  Xt 
by 

C1)  U(t,  X,  Xt)  =  V  X=i>n(t,  X)£n, 

„  ynl 

where  £nJs  are  Wick's  polynomials  (certain  products  of  Hermite  polynomials)  which  form  a  complete 
orthogonal  basis  of  the  probability  functional  space,  and  ^n(i,  x)'s  are  their  deterministic  coefficients 
which  are  also  computed  by  the  inner  product  of  u  and  £n. 

The  most  attractive  point  in  Wiener-Chaos  expansions  is  that  the  randomness  of  a  stochastic 
process  has  been  separated  from  its  deterministic  features.  Therefore,  we  are  able  to  derive  equa¬ 
tions  for  those  deterministic  coefficients  and  solve  them  using  the  standard  methods.  Once  the 
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deterministic  coefficients  ^n(£,  x)Js  are  available,  it  is  easy  to  recover  the  statistics  of  the  solutions 
from  these  coefficients.  For  example,  the  mean  of  the  solutions  is  simple  the  first  coefficient  $)(£,  x) 
in  the  expansion,  and  the  variance  is  computed  as: 

(2)  E(v2)  =  5>a|2. 

Of 

It  has  been  proved  [11]  that  finite  sum  approximations  of  Wiener-Chaos  expansions  is  accurate 
with  an  error  of  the  order  0{e<±tN^1 /{N+l) !),  where  N  is  the  number  of  terms  in  the  Wiener-Chaos 
expansion.  This  suggests  that  if  the  Wiener-Chaos  expansion  is  applied  for  a  short  time  interval 
At,  the  approximation  can  achieve  high  accuracy.  Based  on  this  observation,  we  have  designed  a 
linear  approximation  scheme  (with  iV"  =  0,  which  means  only  one  linear  term  being  used  in  the 
Wiener- Chaos  expansion  at  every  time  interval  with  length  At). 

It  has  been  verified  numerically  that  when  the  variance  of  the  force  term  is  not  large,  espe¬ 
cially  when  the  fluctuation  ||  (vf)2  —  E((vf)2)  ||  is  small,  this  linear  scheme  can  give  very  accurate 
approximations  to  the  statistical  property.  However,  this  linear  scheme  becomes  insufficient  when 
||  (d)2  —  £((02)[|  is  not  small.  In  this  case,  it  is  necessary  to  include  more  terms  (higher  order 
Hermite  polynomials)  in  the  Wiener-Chaos  expansion  at  every  time  interval  to  approximate  the 
solution  accurately.  However,  if  we  use  the  short  time  interval  expansion  as  we  did  for  the  linear 
scheme,  the  number  of  Wiener-Chaos  coefficients  increases  drastically  since  the  nonlinearity  of  the 
equation  brings  all  higher  order  coefficients,  including  the  cross  term  coefficients  between  different 
time  intervals.  This  could  be  computationally  very  expensive. 

To  alleviate  this  difficulty,  we  propose  the  following  remedy.  On  one  hand,  we  would  like  to 
limit  the  number  of  the  Wiener-Chaos  coefficients  by  using  higher  order  Wiener-Chaos  expansion. 
On  the  other  hand,  we  would  like  to  remain  the  accuracy  and  flexibility  of  short  time  interval 
approximations.  To  this  end,  we  divide  time  interval  [0 ,£]  into  two  subinterval  [0,£-5]  and  [t-5,t], 
At  1.  The  first  subinterval  corresponds  to  the  earlier  history  of  the  Brownian  motion  which 

usually  has  small  impact  on  the  solution  at  the  current  time.  The  second  subinterval  corresponds  to 
the  short  time  contribution  to  the  solution,  which  is  not  very  smooth  and  requires  high  resolution. 
In  both  subintervals,  we  construct  Wiener-Chaos  expansions  separately.  Since  there  are  only  two 
subintervals,  we  can  afford  to  compute  the  cross  terms  in  the  Wiener-Chaos  expansions  while 
retaining  some  higher  order  Wiener-Chaos  expansions.  Furthermore,  the  history  term  is  relatively 
small  and  smooth.  Only  a  small  number  of  Wiener-Chaos  coefficients  are  needed  to  obtain  an 
accurate  computation.  Our  preliminary  results  show  that  this  strategy  is  both  efficient  and  accurate. 

(T).  An  Efficient  Domain  Decomposition  Preconditioner  for  Multiscale  Elliptic  Prob¬ 
lems  with  High  Aspect  Ratios 

Many  problems  of  fundamental  and  practical  importance  have  multiple-scale  solutions.  Typi¬ 
cal  examples  include  transport  of  flows  in  strongly  heterogeneous  media  and  heat  conduction  in 
composite  materials.  When  applying  conventional  FEM  domain  decomposition  methods  to  these 
problems  using  linear  or  polynomial  interpolations,  the  convergence  rate  deteriorates  because  the 
coarse  grid  solver  does  not  account  for  fine  scale  heterogeneous  features.  To  attain  a  satisfac¬ 
tory  convergence  rate  it  is  therefore  important  to  construct  a  coarse  grid  solver  which  reflects  the 
small  scale  structures.  Such  a  solver  has  been  developed  by  Hou  et  al,  [6,  7]  who  introduced  the 
Multiscale  Finite  Element  Method  (MsFEM).  The  basic  idea  behind  the  MsFEM  is  to  construct 
base  functions  which  are  adaptive  to  the  local  property  of  the  differential  operator  and  contain  the 
important  subgrid  information. 

An  important  property  with  the  MsFEM  solver  is  that  the  coarse  approximation  space  is  ”  gen¬ 
eralized”  discrete  harmonic  with  respect  to  the  physical  elliptic  operator  that  contains  small  scale 
coefficients.  From  a  theoretical  point  of  view,  this  property  implies  that  the  MsFEM  solver  is  in 
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a  sense  optimal  within  a  general  class  of  coarse  solvers.  Furthermore  it  allows  us  to  interpret  the 
MsFEM  solver  as  a  natural  extension  of  eg.  linear  finite  elements  to  problems  which  include  highly 
oscillatory  coefficients.  In  particular,  if  the  physical  domain  is  partitioned  into  triangles  or  tetrahe¬ 
drons  and  the  elliptic  coefficients  are  quasi-homogeneous ,  i.e.  constant  on  each  coarse  grid  element, 
then  the  corresponding  MsFEM  using  linear  boundary  conditions  to  construct  the  multiscale  base 
functions  simply  reduces  to  standard  linear  finite  elements. 

In  [1]  we  develop,  analyze  and  test  a  nonoverlapping  domain  decomposition  preconditioner  with 
the  MsFEM  solver  as  the  coarse  grid  solver.  The  proposed  preconditioner  falls  into  the  category  of 
Schwarz  methods,  and  the  main  steps  in  our  analysis  is  based  on  the  general  abstract  framework 
for  the  analysis  of  Schwarz  methods.  We  first  demonstrate  that  the  MsFEM  induces  an  ideal 
nonoverlapping  domain  decomposition  preconditioner  in  ID  which  converges  in  one  iteration.  In 
2D  and  3D  the  ability  to  select  proper  boundary  conditions  for  the  multiscale  base  functions 
will  be  important  to  achieve  a  fast  convergence  rate.  In  this  paper,  we  only  consider  the  case 
when  the  MsFEM  solver  is  the  ” oscillatory  extension”  of  linear  finite  elements.  Our  analysis 
shows  that  the  proposed  preconditioner  gives  the  same  rate  of  convergence  for  general  elliptic 
problems  with  oscillatory  coefficients  as  standard  nonoverlapping  domain  decomposition  methods 
using  conventional  coarse  solvers  achieve  for  elliptic  problems  with  quasi-homogeneous  coefficients. 
Moreover,  we  show  that  the  condition  number  has  only  a  relatively  weak  dependence  on  the  local 
oscillatory  coefficients  and  their  aspect  ratio. 

We  perform  a  series  of  numerical  experiments  to  test  the  performance  of  our  preconditioner 
for  elliptic  partial  differential  equations  in  two  dimensions  with  highly  oscillatory  coefficients.  We 
choose  the  coarse  grid  and  the  boundary  conditions  for  the  multiscale  base  functions  so  that  the 
MsFEM  solver  is  the  ”  oscillatory  extension”  of  bilinear  finite  elements.  The  elliptic  coefficient 
function  is  chosen  to  be  the  product  of  a  quasi-homogeneous  coefficient  function  and  a  periodic 
oscillatory  coefficient  function.  We  demonstrate  that  the  MsFEM  induced  preconditioner  proposed 
in  this  paper  shows  a  logarithmic  dependence  on  the  mesh  ratio  H/h  and  is  almost  insensitive  to  the 
local  aspect  ratios  (for  aspect  ratios  as  high  as  1010).  This  confirms  that  the  rate  of  convergence 
of  our  preconditioner  for  general  elliptic  problems  with  oscillatory  coefficients  is  essentially  the 
same  as  standard  nonoverlapping  domain  decomposition  methods  using  conventional  coarse  solvers 
achieve  for  elliptic  problems  with  quasi-homogeneous  coefficients.  We  compare  our  preconditioner 
with  the  preconditioners  obtained  by  replacing  the  MsFEM  solver  with  the  linear  and  bilinear 
finite  element  solvers.  The  convergence  behavior  for  these  preconditioners  deteriorates  rapidly  if 
the  aspect  ratio  within  the  coarse  grid  elements  blows  up.  As  we  are  not  aware  of  any  other  coarse 
solvers  which  successfully  handles  high  aspect  ratios,  the  linear  and  bilinear  finite  element  solvers 
serve  the  purpose  of  illustrating  the  need  for  coarse  solvers  which  can  handle  highly  oscillatory 
coefficients. 

(8).  Singularity  formation  in  3D  vortex  sheets 

One  of  the  classical  examples  of  hydrodynamic  instability  occurs  when  two  fluids  are  separated 
by  a  free  surface  across  which  the  tangential  velocity  has  a  jump  discontinuity.  This  is  called  Kelvin- 
Helmholtz  instability.  Kelvin-Helmholtz  instability  is  a  fundamental  instability  of  incompressible 
fluid  flow  at  high  Reynolds  number.  The  idealization  of  a  shear  layered  flow  as  a  vortex  sheet 
separating  two  regions  of  potential  flow  has  often  been  used  as  a  model  to  study  mixing  properties, 
boundary  layers  and  coherent  structures  of  fluids. 

In  collaboration  with  my  former  Ph.D.  studnet,  Gang  Hu,  [8],  we  have  studied  singularity  forma¬ 
tion  of  3-D  vortex  sheets  using  a  new  approach.  First,  we  derive  a  leading  order  approximation  to  the 
boundary  integral  equation  governing  the  3-D  vortex  sheet.  This  leading  order  equation  captures 
the  most  singular  contribution  of  the  integral  equation.  Moreover,  after  applying  a  transformation 
to  the  physical  variables,  we  found  that  this  leading  order  3-D  vortex  sheet  equation  de-generates 
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into  a  two-dimensional  vortex  sheet  equation  in  the  direction  of  the  tangential  velocity  jump.  This 
rather  surprising  result  confirms  that  the  tangential  velocity  jump  is  the  physical  driving  force  of 
the  vortex  sheet  singularities.  We  show  that  the  singularity  type  of  the  three-dimensional  problem 
is  similar  to  that  of  the  two-dimensional  problem. 

In  the  March  Issue  of  SIAM  News  (2002) ,  there  is  a  long  feature  article  describing  Yizhao’s  work 
with  his  PhJD,  student  in  studying  singularity  formation  of  three  dimensional  vortex  sheets.  The 
article  concludes  that  “the  new  results  of  Hou  and  his  colleagues  have  brought  The  ever  distant 
goal5  of  understanding  turbulence  and  hydrodynamic  instability  at  least  a  few  steps  closer” , 
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